On 



The b Distribution of the Lya Forest: Probing Cosmology and the 

Intergalactic Medium 

Greg L. Bryan 1 
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139 

and 

Marie E. Machacek 



(— | . Department of Physics, Northeastern University, Boston, MA 02115 

oo 

(N 

ABSTRACT 

> 

q^ ■ We investigate a method to determine the temperature-density relation of the inter- 

galactic medium (IGM) at z ~ 2 — 4 using quasar absorption line systems. Using a simple 
model combined with numerical simulations we show that there is a lower cutoff in the dis- 
tribution of column density (Nm) and line width (b parameter). The location of this cutoff 
can be used to determine the temperature-density relation (under certain conditions). We 
describe and test an algorithm to do this. The method works as long as the amplitude of 
fluctuations on these scales (~ 100 kpc) is sufficiently large. Models with less power can 
mimic higher temperatures. A preliminary application is made to data from two quasar 
lines-of-sight, and we determine an upper limit to the temperature of the IGM. Finally, we 
examine the full distribution of 6-parameters and show that this is completely specified by 
just two parameters: the temperature of the gas and the amplitude of the power spectrum. 
Using the temperature upper limit measured with the Nm — b cutoff method, we derive an 

^ ■ upper limit to the amplitude of the power spectrum. 

Subject headings: cosmology: theory, intergalactic medium, quasars: absorption lines 



1 Hubble Fellow 



Introduction 



It has become clear that observations of absorption 



lincfs 



in the spectra oi high-redshiit quasars can give 
us 'valuable information about the nature and distri- 
butio n of the intergalactic medium. Early theoretical 
work (Doroshkcvich fc Shandarin 1977 ; Rees 198q; Bond, 



some recent suggestions on other ways to increase the 
temperature of the IGM. The first is a suggestion of 
Compton heating from a hard X-ray background ( Madau 



& Efstathiou 199£). The second stems from the observa- 



Szalay fc Silk 1988; McGill 199C; Bi, Bonier fc Chu 



1992), supplemented with numerical simulations (Cen et 



tion that the commonly adopted optically-thin limit for 
photo-ionization heating (particularly for helium) may re- 
sult in a substantial underestimate of the gas temperature 
( |Abcl fc Hachnclt 19~99| ). A third, which we will not ex- 
amine in detail in this paper, is provided by photoclcc- 



al. 1994|; |Pctit,jcan, Miiket fc Katcsj; |Zhang ct al. 1995|; tric h eatins from dust grains QNath, Sethi fc Shchckinov 



Hci[nquist ct al. lijjjb j showed convincingly that absorp- 
tion lines at z ~ 6 with column densities less than about 
10 16 cm -2 arise primarily from a network of relatively 
low density filaments and sheets that naturally form out 
of hierarchical primordial perturbations. 

Having established the link between cosmology and the 
Lya forest, subsequent work has focused on two related 
areas: improving our understanding of the physical condi- 
tions of the IGM (at these redshifts), and using the forest 
to constrain cosmological parameters. This includes us- 



ing the power spectrum of the flux distribution (Croft ct 
al. 1998|; |Croft, Hu fc Dave 1999]), the slope of the col- 



umn density distribution (Hui, Gnedin & Zhang 1997; 
Gnedin 1998; Machacek et al. 1999), and inverting the 
flux-density relation (Nusscr fc Hachnclt 1999|). 



Although early simulations seemed to show that all 
models were in agreement with observations, recently it 
was shown (Theuns et al. 1998; Bryan et al. 1999), that 
the width of the absorbers, commonly quantified by the b 
parameter of a Voigt-profile, had been over-predicted in 
most previous work. This left a discrepancy between the 
canonical model and the observations. 

A number of ways to resolve this have been suggested, 
mostly revolving around ways to increase the tempera- 
ture of the gas, and hence the width of the lines. The low 
density gas in the IGM is very close to photoionization 
equilibrium with a background radiation field, usually as- 
sumed to be from quasars. Its temperature is determined 
by a competition b etween adiabatic co oling and photo- 
ionization heating ( Hui fc Gnedin 1997 ). An increase in 
the density will result in more photo-ionization heating 
and hence higher temperatures. In Theuns et al. (1999), 
it was shown that increasing flt>, the ratio of the baryon 
density to the critical density, could widen the lines. How- 
ever, even after doubling the baryon density to the edge 
of that permitted by primordial nucleosynthesis, they still 
found some disagreement. Along similar lines, delaying 
helium rcionization to z ~ 3 — 4 (Haehnelt & Stcinmctz 



1998) can provide a small boost in the temperature. 
In part driven by this discrepancy, there have been 



1999 ). Each of these could, in principle, provide the factor 



of two increase in the temperature required. 

However, since the width of the lines is not just due to 
temperature but also comes from the velocity structure 
(both peculiar and Hubble velocities) along the line-of- 
sight, it seems likely that other parameters also play a 
role. In an elegant paper based on linear perturbation 
theory, Hui & Rutledge (1998) argued that the width 
should depend inversely on the amplitude of the primor- 
dial density fluctuations. 

In this paper we show that there exists a way to indi- 
rectly measure the temperature of the IGM. The method 
is based on a lower cutoff in the Njh — b distribution (first 
noted by Zhang et al. 1997). The position and slope of 
this line is a reflection of the density-temperature relation 
of the IGM. A simple model for this is presented in sec- 
tion and extensive tests using numerical simulations are 
described in section || We develop a simple but robust 
statistic to find the location and slope of the cutoff in the 
Am — b plane. 

However, we also demonstrate that this method breaks 
down if the amplitude of the density fluctuations is too 
low. Again, we present a simple explanation for why 
this occurs and show directly with simulations that it can 
mimic the effect of higher-temperature gas. This means 
that the density-temperature relation derived in this way 
must be treated as an upper limit (until the power spec- 
trum can be fixed by other means). 

Switching from the cutoff in the Ahi — b distribution to 
the full distribution of b parameters, we show in section ra 
that the entire distribution is controlled by the same two 
parameters described above: the temperature of the gas 
and the amplitude of the primordial fluctuations. In fact, 
in this case these two variables are completely degenerate 
and form a single parameter. 

In the last section (section |6|), we apply our tests to 
previously published observations of two quasars, and de- 
rive a temperature density relation. 



2. Theory 

First, we quickly review the calculation of the line pro- 
file; more complete discussions can be found els ewhere 
flHui, Gncdin fc Zhang 1997| ; |Zhang et al. 199^ ). The 



optical depth at a given (observed) frequency i>o can be 
calculated with: 



r M 



n m a a 



dx 

T+z' 



(i) 



where x is the comoving radial coordinate along the line- 
of-sight, and nui is the neutral hydrogen density at this 
point (with redshift z). The Lya cross-section a a is a 
function of the frequency of the photon with respect to 
the rest frame of the gas at position x: 



v = v Q (l + z )(l + ^ 
\ c 



(2) 



Here, v pec is the peculiar velocity and z is the redshift 
due to the Hubble expansion only. This can be rewritten 
in terms of the velocity: 



H 

1 + z 



(x — xo) + v pec (x), 



(3) 



where we are expanding around the point Xq, and H is 
the Hubble "constant" at this redshift. The expression is 
valid as long as u/c is much smaller than 1. In this case, 
the optical depth can be written as: 



T (uo) = Yl 



n-Hi 



1 + z 



du 

dx 



a a du. 



(4) 



The summation sign arises because equation (|3|) can be 
multi- valued. The cross-section, assuming that Doppler 
broadening dominates over natural or collisional line- 
broadening (accurate for column densities less than 10 17 
cm -2 ), is given by: 



c Q ,o 



60F 



-(u-uoY/b 2 



(5) 



Where we have used the following standard definitions: 
b = (2kbT /m v ) 1 ' 2 , kb is the Boltzmann constant, m p is 
the proton mass and <7 Qj o is the Lya line-center cross- 
section. 

2.1. Measuring the temperature of the IGM 

As outlined in the introduction, a central question is 
how to determine the temperature-density relation of the 
gas. We also loosely refer to this as the equation-of-state. 
It is quantified as: 



T = T (1 



■5 h y-\ 



Here, (1 + 5b) — pb/(£lbf>), where p is the mean density 
(both baryonic and dark). For gas which is primarily 
heated by UV photoionization, 7 is expected to vary from 
1 immediately after reionization, to a limiting value of 
about 1.5 (Hui & Gnedin 1997). Similarly, T is expected 
to evolve as To ~ (1 + z) lJ . This information is encoded 
in the absorption lines, and here we describe a way to 
indirectly measure, or at least constrain, the equation of 
state. 

The width of a given line is the result of a convolu- 
tion involving both the temperature and velocity distri- 
butions of the gas along the line of sight. However, for 
low column-density lines (< 10 15 cm -2 ) both the density 
and temperature are slowly varying functions of position 
(Bryan et al. 199E). Therefore, it makes sense to distin- 
guish two sources of the total line-width: 

• &t, the thermal Doppler broadening, which is a 
measure of the optical-depth weighted temperature 
of the gas, and 

• b ve i, the bulk velocity broadening, which in turn 
has two components: the peculiar velocity and the 
Hubble velocity. 



(6) 



The total width comes from adding these two in quadra- 
ture: btotal = (&| + &L) 1/2 - 

Our method is based on two assumptions. First, that 
the column density of a line is proportional to the density 
of the gas, so a measurement of iVui can be converted to 
(1 + Sb). The second assumption is that there are at least 
some lines (at a given TVhi) for which b ve i is significantly 
smaller than br so that btotal ~ &r- If this is true, then 
there will be a minimum in the TVhi — b distribution given 
by: 

b min ~b T = {2k b T/m p ) 1 / 2 (7) 

= (2k b To(l + S b r- 1 /m p ) 1 / 2 : 

where 1 + Sb is the mean overdensity of a line with column 
density iVui. In fact, this argument was first suggested 
by Zhang et al. (1997). 

To make this a little more concrete, we can gener- 
ate a toy model for this minimum based on a number 
of assumptions: (1) photoionization equilibrium holds, 
(2) when computing column densities, peculiar velocities 
can be ignored, (3) all systems have the same comoving 
length, and (4) at any column density, there exist some 
absorbers with b ve i = 0. None of these conditions hold ex- 
actly, however numerical simulations show that — for at 



least some models — they are not unreasonable approxi- 
mations ( |Bryan ct al. 1999 ). In fact, as we will show, it is 
the last assumption that will present the most difficulties. 
From the first assumption, it is relatively straightfor- 
ward (e.g. Zhang et al. 1998) to show that the neutral 
hydrogen density is given by: 



n H i = 1.2xKr ib (l + <5 6 ) 



0.02 



2 (i+fl! : 

ThL-12 



density fluctuations on the scales giving rise to the for- 
est. The most convincing demonstration of this comes 
from numerical simulation; however, a simple plausibility 
argument can be made as follows, fj 

We begin with a sinusoidal perturbation of comov- 
ing wavenumber k which perturbs a fluid element's La- 
grangian position q: 



-°- 7 cm- 3 , 



(8) 
where T4 — T/10 A K and 1^1,-12 is the hydrogen pho- 



D + Asin(kq)/k, 



(13) 



toionization rate in units of 10 



-12 



The next two 



assumptions provide a relation between this and the col- 
umn density, which we will take to be: Njjj — Inm where 
I = 125 (4/(1 + z)) kpc, so that: 



N, 



HI 



1.9 x 10 lci F(l-|-(5 b ) 2 rf u 'cni 



(9) 



For notational ease, we have taken the cosmological and 
photo-ionization factors into F which is defined as 

f=[^-\ (±±£) r-, 1 10 . no) 



where A is the initial amplitude and D + describes the 
evolution of a growing mode for the given cosmology 
(D + oc (1 + z) _1 in an Einstein-deSitter universe). This 
equation usually begins a discussion of the Zel'dovich ap- 
proximation but here we will need to assume that A is 
small compared to unity. 

The peculiar velocity of the fluid element is given by 



r 



pec 



= ax — afD + Asm(kq)/k. 



(14) 



0.02 



12- 



In fact, the exact value of I has been selected to give a 
good fit to the data described below, however it is quite 
compatible with the width of filaments seen in simula- 
tions. 

Finally, we use the fourth assumption along with equa- 
tion (pi) to derive an expression for the minimum column 
density as a function of temperature: 

N HI . mm = l.gxio 13 ^ 2 /^- 1 )- - 7 ^ 2 ^ 7 "^^- 2 , (11) 

where T 0> 4 is T /10 4 K. Using equation (R) this can be re- 
cast entirely in terms of observable quantities (NHi,min, b m in 
and parameters of the equation of state (To, 7): 



We adopt the common notation / = aD + /aD + , where 
a = (1 + z)^ 1 is the scale factor ( Peebles 1993 ). For 
gas in photo-ionization equilibrium, the neutral hydrogen 
density is related to the gas density by nni ~ p 1 ' 7 ■ The 
density produced by this perturbation is given by p ~ 



(1 



D + Acos(kq)) so, 



Tim oc (1 + D + Acos(kq)) oc (1 - 1.7D + Acos(kq)) . 

. (15) 
For brevity, we have dropped the coefficients to this ex- 
pression since they contribute only to the overall normal- 
ization of the optical depth, not the structure of the line. 
This is the density in physical space. In order to 
' compute the redshift-space density and hence the optical 
depth distribution via equation (0) , we need the Jacobian, 



13km/s 



1.9 x 10 13 cm" 2 F 



T, 



1/(2.7-0.77) 



0,4 



(12) 

This expression shows that the minimum in the Nai — b 
distribution should take the form of a power law (since 
the equation of state is assumed to be a power law). It 
gives a way to determine the parameters of the density- 
temperature relation in equation (o) from a measurement 
of the intercept and slope of the minimum line in the 
^Vhi — b plane. 

2.2. The cosmology-6 connection 

What else, besides temperature, influences the width 
of the absorbers? The most significant cosmological pa- 
rameter turns out to be the amplitude of the primordial 



dv 
dx 



-(l + fD + Aco S (kq)). (16) 

a 



Using the previous two expressions, the full expression for 
the optical depth distribution is given by: 

rex Ml- (1.7 - f)D+Acos(ku/a))a a du, (17) 

where we have employed equation (0) (to first order in A) 
to write an expression only in terms of u. 



2 Another calculation along these lines, but for a random Gaussian 
field instead of a single perturbation, can be found in Hui & Rut- 
ledge (f998), where they also derived the expected shape of the 
distribution of b parameters. 



While this is helpful, there is still a convolution with 
the Doppler width to contend with. The general expres- 
sion is quite complicated; it is more helpful to recognize 
that the result of the convolution will be quite close to a 
Gaussian with width {b^+b 2 ^ 1 / 2 (see Bryan et al. (1999) 
for an explicit demonstration of this). The first term is the 
Doppler broadening contribution while the second term 
comes from the structure of the line, as given in equa- 
tion (pjj). Under the assumption that these two terms are 
independent, we can ignore the thermal broadening part 
and focus simply on the velocity part. However, the si- 
nusoidal perturbation is not Gaussian, so determining a b 
parameter from this is not trivial. Of course, this is quite 
true in the real forest, where line profiles are often not 
well described by Voigt-profiles. We make the correspon- 
dence by matching the shapes of the cosine and Gaus- 
sians profiles near the peak of the line, where the largest 
contribution to the total optical depth occurs. This is 
done by expanding both the cosine term of equation ([17]) 
and a Gaussian exp(— u 2 /b 2 ) and equating the first non- 
constant term, which is proportional to u 2 . Another way 
to do this would be to start with a Gaussian perturbation 
instead of the sinusoidal one in equation (Uq). In either 
case, we find that: 



h 2 

°vel 



2H 2 



(1.7- f) AD +k 2 (l + z) 2 ' 



(18) 



It is interesting to examine this expression in more 
detail. For the redshifts under consideration here, it 
is useful to approximate the Hubble velocity as H = 
Hq^Iq (1 + z) 3 ' 2 , which is accurate as long as Qq, the 
ratio of the total matter density to that required to close 
the universe, is not too low (Peebles 1993). Similarly, 
f rs J7 06 sw 1, since J7 is close to one at this redshifts 
(again, as long as Q,q is not too low). By the same rea- 
soning, the growth factor D + « (1 + z)~ x . 

The wavenumber of the perturbation k is also a factor 
in this expression. While there will be a range of wave- 
lengths, it seems very reasonable to associate this with 
the Jeans wavenumber or at least some fixed fraction of 



it (e.g. Hui, Gncdin fc Zhang 1997; Gncdin fc Hui 1998 



kj — 



12ira 2 Gpp 
5fchTo 



(19) 



In this expression, \i is the mean mass per particle and 
is about 0.6m p for ionized gas. The average density p = 

3H$Q (1 + z) 3 /8ttG so kj oc H^ 1 ^ 1/2 . (In writing this 
we have suppressed a factor of (1 + z) _1 Tq 5 which is 



likely to be quite small since ro~(l + z) 17 ). Therefore, 

k oc Hq S7 (1 + z) since it is a comoving wavenumber. 

Surprisingly, if we make these assumptions then Hq, Qq 

and z all cancel, leaving the remarkably simple expression 



el OC 



.4 



-0.5 



(20) 



For a more realistic model with a spectrum of pertur- 
bations -P(fc), the amplitude A will be proportional to 
the amount of power on the scales of interest (i.e. A 2 oc 
P(kj)). For a given spectral shape b ex <r 8 . In fact, 
approximately this scaling was found in Machacek et al. 
(1999). It fails when the perturbations become too large 
or when the thermal width dominates for a majority of 
lines; however, it works surprisingly well, even in the 
trans-linear regime. It also explains why the other cosmo- 
logical parameters, in particular f^o and the Hubble con- 
stant, have so little effect on the b parameter distribution. 
The lack of a redshift dependence in equation ( |20| ) also 
helps to explain why the median distributions (both simu- 
lated and observed) seem to be so constant with redshift, 
while a thermally dominated distribution would scale as 
b ~ T 1 ' 2 - (1 + z) - 8 . We note that this result dif- 
fers slightly from that derived in Hui & Rutledge (1999), 
who found b oc (D+erg) -1 / 2 because they assumed that 
the smoothing scale would be constant in redshift space 
rather than comoving space. 

3. Simulations 

In order to examine these effects in more detail, we 
have performed numerical simulations of a range of mod- 
els with various heating rates and hence various equations 
of state. We use a grid-based method based on the piece- 
wise parabolic algorithm to model the gas and a particle- 
mesh code for the dark matter and gravity. We follow 
the abundances of six species: HI, HII, Hel, Hell, Helll 
and e _ by solving the non-equilibrium evolution equa- 
tions. The simulation method is described in more detail 
elsewhere (|Bryan et al. 1995; Anninos et al. 1997). 



Table [P lists the simulations that we analyze in this 
paper. The first column gives the cosmological model — 
most are a form of the currently popular cosmological con- 
stant dominated (LCDM) with fi = 0.4, fl A = A/3-ff, 2 = 
0.6, fif, = 0.05, and h — 0.65 (the Hubble constant in 
units of 100 km/s/Mpc). We also run one other model in 
order to demonstrate that the cosmological dependence 
is well-understood. This is a flat model (SCDM) with 
ilo = 1, fib = 0.08 and h = 0.5. The second column 
shows the selected normalization of the power spectrum 



by giving erg, the linearly extrapolated rms density fluctu- 
ations in a tophat sphere of 8ft. -1 Mpc. The next column 



giv ;s a measure of small scale fluctuations suggested by 
Gnedin (1998): 



tional volume. These spectra are then analyzed with 
an automated Voigt-profile fitting routine (Zhang et al 
1997). This algorithm does not include a number of ob- 



'34 



P{k,z = 3)e 



-2k 2 /k 2 3 



k 2 dk 

: 

2tt 2 



(21) 



servational effects such as noise and so is somewhat ideal- 
ized; however a comparison with a more realistic method 
indicated that for high signal-to-noise spectra the differ- 



where P(k, z = 3) is the power spectrum at z = 3 and 



ences are not large (Bryan et al. 1999); we will discuss 
this point in more detail below. 



S:M 



,1/2 



34£1 / h Mpc 1 . All power spectra in this paper 4. The density-temperature relation 



come from the analytic fits of Eisenstein & Hu (1999). 

The fourth column indicates L the size of the simula- 
tion vo lume, in Mpc. As we demonstrated in a previous 
paper ( Bryan et al. 1999 ), there is some dependence of 
the 6-parameter distribution on the box size, since fluctu- 
ations with wavelengths larger than L are not included. 
Our canonical box size of 4.8 Mpc is sufficient for a reason- 
able prediction, however convergence requires 9.6 Mpc, so 
we perform one simulation with this larger size. All runs 
use a grid of 128 3 cells (except the 9.6 Mpc box which uses 
256 3 ); this provides the minimum resolution required to 
accurately resolve the line profiles. 

The radiation field is assumed to be spatially constant 
with the form given by Haardt & Madau (1996), which 
assumes that the ionizing photons come from the observed 
quasar distribution. However, we modify the Hell photo- 
heating rate in order to account for the neglected radiative 
transfer effects as discussed in Abel & Haehnelt (1999). 
Although this is not realistic in detail it does produce 
the desired result of heating the IGM. The fifth column 
of Table [j] indicates the factor by which this is increased 
relative to the original Haardt & Madau heating rates. 
Abel & Haehnelt suggested this factor should be ~ 2 — 4. 

The next column indicates if the simulation includes 
Compton heating due to a hard X-ray background. We 
use the heating rate as computed by Madau & Efstathiou 
(1999) who assumed that the energy density evolved as 
U x {z) = Ux{0){1 + z) 4 cxp(-z 2 /z2) and included the 
Klein- Nishina relativistic corrections to the cross-section, 
resulting in a heating rate that scales approximately as 
(1 + z) 13 / 3 . We adopt z c = 5. 

The next four columns indicate the equation of state 
parameters as given in equation (0), at z = 2.7. The 
first set (without primes) come from fitting the Nm — 
b minimum and are the observational estimates, while 
the second set (with primes) come from directly fitting 
10,000 randomly selected cells in the simulation. The 
last column is the median value of the b distribution. 

The analysis is carried out by generating artificial spec- 
tra along random lines-of-sight through the computa- 



4.1. Testing with simulations 



In section 2.1 



we motivated why there should be a 
minimum in the Nm-b bivariate distribution. In Fig- 
ure EL we show this distribution for our LCDM run with 
<j$ = 1-0 and the usual Haardt & Madau (1996) photo- 
heating rates for three redshifts: z = 4, 3 and 2. There is 
a sharp cutoff at low column densities and large b values 
(i.e in the upper-left corner of each frame) which is due 
solely to our criterion for identifying lines, namely that 
the optical depth at the line center be larger than 0.05. 
More interestingly, there is another fairly sharp cutoff at 
low b which is the subject of this paper. 

The sharp edge that defines the cutoff is fairly obvious 
to the human eye a nd has been previously no ticed in both 
observations (e.g. Kirkman fc Tytler 1997 ) and simula- 



tions ( Zhang et al. 1997 ). In order to be more quantita- 
tive about the position of the cutoff, we take as our in- 
spiration edge-detection techniques from machine-vision 
research. An edge — in one dimension — is defined as a 
zero in the second derivative of the intensity (since this 
is an extrema in the first derivative, the rationale is obvi- 
ous). The application of this idea is quite straightforward. 
First, we sort the lines by column density and divide 
them into groups of size 30-50 (each group is equivalent 
to a scan line in an image). The smoothed density of lines 
is then computed as a function of b with a weighted sum 
over all lines in each group: 



Pb(b) 



cxp(-(6 J -6) 2 /2a 6 2 ), 



(22) 



where Ob = 3 km/s is the smoothing constant. The 
method is not sensitive to small changes in either this 
parameter or to the number of points in the group (more 
lines per group mean less noise but lower resolution along 
the Nm direction). We can compute derivatives of pb very 
easily, so we simply define, for each group with average 
column density Nhi minj the edge to be at 6 m ; n such that 



(Ppb 
db 2 



\pmin) U. 



(23) 



Table 1 
Simulations analyzed in this paper 



model 


0"8 


0"34 


L (Mpc) 


rHcIl/rHcII.HM 


X-ray heating? 


To, 4 


7 


-^0,4 


V 


Vcd (km/s) 


LCDM 


1.0 


1.93 


4.8 


1.0 


no 


1.01 


1.40 


1.07 


1.44 


20.1 


LCDM 


1.0 


1.93 


4.8 


2.0 


no 


1.30 


1.39 


1.35 


1.43 


23.0 


LCDM 


1.0 


1.93 


4.8 


4.0 


no 


1.80 


1.37 


1.84 


1.40 


27.1 


LCDM 


1.0 


1.93 


4.8 


1.0 


yes 


1.18 


1.34 


1.28 


1.26 


21.5 


LCDM 


0.8 


1.54 


4.8 


2.0 


no 


1.52 


1.34 


1.36 


1.42 


24.6 


LCDM 


0.8 


1.54 


1.8 


1.8 


no 


1.44 


1.37 


1.39 


1.26 


23.9 


LCDM 


0.8 


1.54 


9.6 


1.8 


no 


1.49 


1.32 


1.42 


1.27 


24.3 


LCDM 


0.6 


1.16 


4.8 


2.0 


no 


2.31 


1.22 


1.44 


1.32 


28.9 


SCDM 


0.55 


1.55 


4.8 


2.0 


no 


1.53 


1.32 


1.29 


1.38 


24.6 


LCDM/iO.5 


0.8 


1.54 


4.8 


2.0 


no 


1.22 


1.40 


1.17 


1.43 


24.3 
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Fig. 1. — The column density- width (Nui—b) distribution 
at three redshifts (z — 4, 3, 2) for our LCDM model with 
as = 1.0 and our standard heating rates. The points are 
about 1200 simulated lines, the filled diamonds trace the 
minimum in the JVhi — b distribution as described in the 
text, while the solid line is a fit to these points. 



For noisy data there are occasionally several zero cross- 
ings — we take the strongest, defined as the one with the 
largest first derivative. In order to get the lower cutoff, 
we insist that the first derivative be position. Software to 
perform this is available from the authors on request. 

One important advantage of this algorithm is that it is 
relatively insensitive to noise which tends to smear out the 
data (but not shift the edge) or outliers in the Nja — b 
distribution (which may be rogue metal lines or simply 
the result of blending) . It also non-parametric in that it 
doesn't assume a form for the iVui — b line. 

This results in a set of points which define the TVhi — b 
minimum and are plotted in Figure |lj as solid diamonds. 
We also show a least-squares power-law fit to the line. 
We used only absorption lines within a range of column 
densities which was selected to include the majority of 
lines but not go below about 2 x 10 12 cm -2 or above a few 
times 10 15 cm~ 2 . The lower limit is slightly below present 
day observational limits and the upper limit marks the 
point were line dynamics becomes more complicated (i.e. 
affected by shocks and, in the real world, star formation). 

Each of these points may be converted into a measure- 
ment of T and 5b via equations (||) and (0) . The results 
are plotted in Figure || as open circles. Also shown is a 
measure of the temperature-density relation in the simu- 
lation by plotting 10000 random cells as dots. The solid 
line comes from converting the power-law fits from Fig- 
ure into measurements of To and 7 via equation (fl3) . 

The match is quite good, although clearly we prefer- 
entially probe the upper part of the T — 8 relation. This 
is because, to be observed, a line must be more dense 
than the surrounding gas, so the measure is insensitive to 
the temperature of gas between filaments and sheets (al- 



though note that we can still probe densities significantly 
lower than the cosmic mean) . On the other end, the max- 
imum overdensity is around 10 because of the maximum 
TVhi limit. If the T — 5 relation is not a strict power-law, 
as at z = 4, then this can result in a substantial under- or 
over-estimate of the temperature at very low or very high 
densities. We also note that there are a small fraction of 
points in the simulation with moderate densities but very 
high temperatures. These points tend to lie near much 
more massive structures and have been enveloped in their 
accretion shock. The minimum b method described here 
is not sensitive to these (rare) points. 

We should remind the reader of two points about the 
normalization of the 5b — TVhi relation, equation (|12|). 
First, the normalization was selected to give a good match 
at z — 3 — changing this value is roughly equivalent to 




Fig. 2. — The temperature-overdensity relation for the 
same canonical LCDM model (as = 1.0) in Figure EL 
shown at three redshifts: z — 4,3 and 2. The dots are 
10000 random cells, while the open circles and solid line 
are derived from the minimum Ahi — b line as described 
in the text. 



shifting the open circles horizontally (in 1 + 5b). Second, 
the parameter F shows there is a degeneracy amoung the 
parameters Qj,, h and Thi such that as long as the value of 
F is unchanged, these parameters can be changed with- 
out affecting the results plotted here. In fact, this is one 
reason we choose to plot 1 + 5b rather than a physical 
density. Although the individual parameters f2&, h and 
Fhi are not well determined, this particular combination 
is, from observations of the Lya forest (see, for example, 



Rauch et al. 1997) 



4.2. Changing the equation of state 

Given our success in measuring the temperature of the 
IGM in the canonical simulation described in the previ- 
ous section, it is interesting to see if we can detect the 
effect of changing the primary heating mechanism as out- 
lined in the introduction. In this section, we show that 
this is possible. We retain the same cosmological model, 
but modify the Hell photo-heating rate as described in 
Abel & Haehnelt (1999). This takes into account radia- 
tive transfer effects during Hell photo-ionization that are 
neglected in these simulations. Since the amplitude of the 
effect is hard to gauge, we either multiply the rate by two 
or four. 

The results are shown in Figure 3, again with b m i n 
points and fit as determined from our edge-detection al- 
gorithm. Here, in order to give a concrete comparison 
with observations and to provide a constant reference 
point, we also plot the observational equivalent of the 
Nui—b minimum as found by Kirkman & Tytler (1997) at 
a mean redshift of 2.7 for a single line-of-sight. Although 
they fit this by eye the result agrees very well with the 
method used here. Our standard Hell photo-heating rate 
produces temperatures which are too low, while the x2 
and x4 simulations are much closer and bracket the re- 
sult, with the x2 case the closest. There is some evidence 
that the slope is in disagreement for all cases; however, 
this does not appear to be strongly significant given our 
level of uncertainty (see below) . We can apply the same 
method used earlier to determine the T—5 relation, which 
is shown in the three left panels of Figure |[ Again, the 
equation of state is quite accurately determined. 

Next, we examine the importance of Compton X-ray 
heating in the bottom panel of Figure |3|, which is again 
our canonical LCDM simulation with the usual Hell 
photo-heating rates. However, now we include hard X- 
ray heating as described earlier. While this does boost 
the temperature somewhat, it is clearly — by itself - 
insufficient to match observations. Since the Compton 
heating rate is independent of density, it tends to flat- 
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Fig. 3. — The column density-width (iVni — b) distribu- 
tion at z — 2.7 for our canonical LCDM model with four 
different heating rates. The top frame shows the usual 
Haardt & Madau (1996) Hell photo-heating rate, while 
the next two panels demonstrate the effect of increasing 
this rate by factors of two and four, respectively. The bot- 
tom panel has the usual Hell rate but includes Compton 
X-ray heating. The filled diamonds describe the mini- 
mum in the AThi — b distribution as described in the text. 
The solid line is a fit to these points, and the dashed line 
is the same for each panel and shows the observational 
determination of the Nm — b minimum from Kirkman & 
Tytler (1997). 



ten the temperature-density relation, shown in the upper 
right panel of Figure [|. However, the effect at z = 2.7 is 
mostly limited to low densities and so is very difficult to 
detect with the minimum TVhi — b method. 

4.3. Changing the power spectrum 

Although we have been successful in measuring the 
equation of state for our canonical LCDM model, we ar- 
gued in section that the amplitude of fluctuations on 
scales of a few hundred kpc is also important in determin- 
ing the distribution of line widths. In this section, we will 
demonstrate that for some models, this effect prevents us 
from accurately measuring the temperature-density rela- 
tion. 

Figure shows the results from two LCDM models 
in which the power has been reduced to a% = 0.8 and 
as = 0.6. This changes the derived minimum Ahi — b 
line. For ag, = 0.8 (the top panel) the match with ob- 
servations is very good, while the lower power simulation 
produces a power-law fit that is too flat. Since the equa- 
tion of state has not changed, application of the iVni — b 
minimum method results in temperatures which are too 
hot, particularly for the lowest-power run. This shown in 
the bottom right panels of Figure [|. 

This happens because one of our key assumptions is vi- 
olated: specifically that there be a substantial number of 
lines for which b ve i be small. For the lower-power models, 
the smallest filaments are almost all in the pre-turnaround 
stage. That is, their peculiar velocities are still smaller 
than the Hubble flow across their width, so that these 
two values cannot cancel. This preferentially affects low 
column density lines because they are the smallest fluc- 
tuations. In fact, lines with a column density of around 
TVhi ~ 10 15 cm" 2 (1 + 5 b ~ 10) faithfully reproduce the 
correct temperature even for the as = 0.6 simulation. 

This demonstrates that the minimum 7Vhi — b method 
suffers from a degeneracy between the gas temperature 
and the amplitude of fluctuations. As long as CT34 <; 1.6, 
the method can be used in a straightforward fashion (we 
use 0-34 since this is much closer to the scale and redshift 
of interest and so is nearly independent of other cosmo- 
logical parameters). Below this value, there is still use- 
ful information to be gained, but the interpretation is 
more complicated. In particular, without other knowl- 
edge about CT34, the value of Tq derived in this way is an 
upper limit and the value of the slope 7 is a lower limit. 

Although we do not show the results here, we have also 
analyzed an LCDM simulation with a lower value of the 
Hubble constant as well as an SCDM model (see Table 
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Fig. 4. — The temperature-overdensity relation at z = 2.7 for six simulations. The symbols are as in Figure M. The 
three left panels are from the canonical LCDM model (with as = 1.0) for the same three values of Hell photo- heating 
rate as in Figure 0. The upper right panel is for the same model but with X-ray Compton heating instead. The two 
lower panels on the right were run with twice the usual Haardt & Madau Hell photo-heating rate but for low-power 
LCDM models, with erg = 0.8 and a$ = 0.6. For comparison, each plot also shows the equation of state from the 
middle left frame as a dashed line. The triangles and dot-dashed line in the middle-left panel show the equation of 
state derived from a different Voigt-profile fitting technique as described in the text. 
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1). The results agree with the trends discussed in this 
section. 

It is important to ask at this point what the uncertain- 
ties are in determining the minimum Nm — b line. The pri- 
mary source of uncertainty is fitting the Voigt-profilcs in 
the first place, as this is both non-linear and non-unique 
(e.g. Kirkman & Tytlcr 1997). In order to gauge the 



magnitude of the possible error, we select one simulation 
(LCDM eg — 1.0 with twice the usual Hell photo- ionizing 
rate) and fit it with the more realistic method AUTOVP 
(Dave et al. 1997), kindly provided by Romeel Dave. This 
method performs a chi-squared minimization to produce 
the line list from the simulated spectrum. Figure ra shows 
the result for a signal-to-noise ratio of 60, along with the 
fit found with the more idealized Voigt-profilc fitting algo- 
rithm. Clearly there is some difference, which also affects 
the derived equation of state, shown in the middle-left 
panel of Figure [|. This amounts to about a 15% difference 
in T at St = 0, and is mostly due to fitting Voigt-profilcs 
to lines which do not follow this profile in detail. 
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Fig. 5. — The column density-width distribution (iVjji — 
b) at z = 2.7 for two LCDM simulations. These both 
use twice the usual Hell heating rate, but decrease the 
amplitude of the initial fluctuation spectrum to rjg = 0.8 
(top) and (78 = 0.6 (bottom). The symbols are as in 
Figure g. 



5. The median of the 6-distribution 

We have so far focused on the low b cutoff, but now 
turn, briefly, to the rest of the distribution. The top panel 
of Figure shows dn/db at z — 2.7 for our three ag = 
1.0 models with varying equations of state (i.e. different 
Hell photo- heating rates) . Only lines in the range Ahi = 
10 13 ' 1 — 10 14 cm -2 are used so as not to be biased by the 
line-selection function. All distributions are normalized 
so that J(dn/db)db = 1. This plot shows that increasing 
the temperature results in a constant shift in log(6). As 
discussed in Bryan et al. (1999), this is not primarily a 
result of thermal Doppler broadening, but comes instead 
from a thickening of the filaments and sheets (in both 
physical and velocity space) due to the influence of the 
increased pressure — gas is driven out of the centers of 
the filaments. 

The middle panel of Figure |7] shows the effect of chang- 
ing the amplitude of the power spectrum while keeping 
the temperature constant. The results appear to match 
the simple scaling derived in section 0. This degeneracy 
between temperature and power can be written in terms 
of the median of this distribution: 

1/2 



b me d(TQ,(T34) = 26.5km/s 



f_a_y 

Viooooxy 



(a 3i y 1/2 (24) 



We use Tq to indicate the temperature at 6b = mea- 
sured directly from the simulations, rather than via the 
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Fig. 6. — The column density- width distribution (iVki— b) 
at z = 2.7 for a LCDM simulation with a% = 1.0 and 
twice the usual Hell photo-heating rate. This is the same 
simulation analyzed in the second panel of Figure pL but 
here we use a different Voigt-profile fitting technique. The 
diamonds and solid line are the derived minimum Nm — b 
line, while the dot-dashed line is the minimum Njh — b fit 
found previously, for comparison. 
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minimum Ahi — b method. In Figure g, we plot this func- 
tion against the measured median b parameters from our 
simulations. Finally, in order to demonstrate that there 
is not much change in the shape, we rescale the b distribu- 
tions in the top two panels of Figure R with the following 
transformation: b — > fb and (dn/db) — > (dn/db)/ f . We 
define / = b med (T^ <r 34 )/b med (13000K, 1.93), where b med 
is given in equation (J24|) . The result is shown in the bot- 
tom panel of this figure. The scaling works very well, in- 
dicating that the shape of the distribution changes little, 
once the median is specified. The biggest differences are 
at small b, where thermal Doppler broadening dominates. 
Our simulated boxes are too small to fully contain all 



the large-scale power. Previously (Bryan et al. 1999) we 
showed that this has little effect on the shape of the b 
distribution but can cause a small shift in the median. 
In order to gauge the size of this effect here, we ran two 
models with ag = 0.8, one with our usual box-length of 
4.8 Mpc and one with twice this size (see Table El). The 
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Fig. 7. — The simulated b distribution function for three 
models with the same power but different temperatures 
(top), the same temperature but differing power (middle) 
and all models scaled as described in the text (bottom). 



both had the same cell size (i.e. the larger box simulation 
had 256 3 cells rather than our more usual 128 3 ). The 
change in the median was only 0.4 km/s (1.6%). This 
could be slightly larger for models with more power, but 
is unlikely to be a significant effect. 

6. A preliminary comparison to observations 

In this section, we make a preliminary comparison to 
observations, using previously published results from the 
quasar HS1946+7658 ( JKirkman fc Tytlcr 1997|) at a mean 
redshift < z >= 2.7 and APM 08279+5255 flEllison ct al" 



1999 ) at < z >— 3.4. Both observations have high signal- 



to-noise ratios, ranging from 15 to 100 for HS1946+7658 
and 30 to 150 for APM 08279+5255. It should be kept in 
mind that the Voigt profile fitting technique used in these 
papers (which is not fully automated) differs somewhat 
from both methods used here. The ATjji — b distributions 
are shown in Figure H along with the derived minimums 
using our edge-detection method, (for APM 08279+5255 
we adopt a minimum column density of 2 x 10 13 cm^ 2 due 
to concerns about line-blending at these high redshifts). 

The minimum Ahi — b edges detected can be converted 
into a measurement of the temperature-density relation, 
which is shown in Figure nOl The values of F required 
to convert column density to 1 + Si, were determined by 
fitting the column density distribution to the simulations 
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Fig. 8. — The median of the b distribution depends almost 
entirely on just two parameters: the temperature of the 
gas (Tq) and the amplitude of the power spectrum (034). 
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(we use the LCDM eg — 0.8 simulation but this is not 
very sensitive to which model we select). The noise in the 
T — 5b relation is larger than for the simulations because 
of the much smaller number of lines (~ 300 as compared 
to ~ 1200). The power law fits are given by To. 4 = 1-65, 
7 = 1.29 (HS1946+7658) and T 0A = 1.52, 7 = 1.54 (APM 
08279+5255). We remind that reader that these determi- 
nations are really upper limits to the temperature rather 
than measurements because of the possible effects of cos- 
mology (i.e. the unknown value of CT34). 

There is an indication from the higher-redshift system 
that the gas is cooler at low density (i.e. a steeper equa- 
tion of state), although clearly this is substantialy uncer- 
tainty. Still, a similar trend of lower b lines at higher red- 
shi ft has b e en previously no ted from different data (Hu 



et fll. 1995| ; [Kim et al. 1997J ), so it is worth considering 



the possibility that (low density) gas is cooler at higher 
redshift. This does not agree with what is expected for 
gas dominated by steady photo-ionization heating and 
adiabatic cooling: Tq ~ (1 + z ) 1 - 7 with a slowly steep- 
ening equation of state slope ( Meiksin fc Madau 1993 
Miralda-Escude & Rees 19941; |Hui & Gnedin 19971; lAbel 



& Haehnelt 1999 ). An alternate heating source, such as 
late Helium reionization, would be required if this result 
proves true. 



100 




Fig. 9. — The column density- width distribution (JVhi — b) 
at for two sets of observed absorption line systems. The 



top is HS 1946+7658 (Kirkmanfc Tytlcr 1997) at < z > 



2.7 and the bottom is from APM 08279+5255 flEllison et 
al. 1999|) at < z >= 3.4. The symbols are as in Figure 



To compare the shape of the b distributions to obser- 
vations, in Figure O we plot dn/db from the same two 
quasar systems previously discussed. We also show the 
distribution from our LCDM er 8 = 0.8 simulation with 
twice the Haardt & Madau Hell photo-heating rate. The 
shape is in reasonable agreement with HS1946+7658, but 
the other system has a large number of very low b lines 
which are not seen in any of the models considered here 
(although a very low temperature model might match). 
Both observed systems also have a more pronounced non- 
Gaussian tail at large b than appears in the simulations. 
Note, however, that the log-log plot accentuates this tail 
when compared to the more usual linear plot. It should 
also be kept in mind that the Voigt-profile fitting method 
used for the observations differs from either employed in 
this paper. Clearly, a more definitive result will require 
identical treatment of data and simulations. The bot- 
tom panel of the same Figure shows how the two different 
Voigt-profile fitting algorithms used in this work compare. 

The observed median b (in the same column density 
range considered earlier) is 27.3 for both systems. Using 
equation (124), this implies that 



(T^/lOOOOif) 172 (034) V2 = 1-03 



(25) 



If we use the upper limits for Tq derived earlier from the 
Ahi — b distribution, then we can get an upper limit on the 
amplitude of the power spectrum: 034 <J 1.52. The uncer- 
tainty is about 20%, due to the uncertainty in the mea- 
sured temperature. This is quite close to the minimum 
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Fig. 10. — The temperature-density relation as derived 
from the two quasars analyzed in Figure M. 
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value of 034 required for a straightforward interpretation 
of the iVfji — b minimum method (1734 ;> 1.6). Interest- 
ingly, this value of 034 (around 1.5-1.6) agrees well with a 
number of determinations of the power spectrum ampli- 
tude using other characteristics of the Lya forest QGncdin 
199$ |Croft et al. 1999| ). For the LCDM model it is also in 
accordance with the normalization from COBE and rich 



clusters of galaxies (e.g. Liddle et al. 1996) 



7. Conclusions 

In this paper, we have investigated a method for de- 
termining the relation between density and temperature 
(loosely denoted the equation of state of the IGM) from 
the distribution of quasar absorption lines in the ./Vhi — b 
plane. Specifically, we look for a sharp minimum line 
in this plane arising from the fact that Doppler thermal 
broadening sets a minimum line width. Because there is a 
tight relation between column density and overdensity, we 
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Fig. 11. — Top: b distributions from two observed quasar 
line systems along with one of the better fitting simulated 
results; bottom: dn/db for the same simulation using two 
different Voigt-pr ofilc fitting algorit hms: the solid line 
shows ZANM97 ( ^hang et al. 1998[ ) which is quit e ide- 
alized, while the dashed line indicates AUTOVP (Dave 



can relate ./Vhi to overdensity and b to temperature. We 
derive a simple model which reproduces this behaviour, 
and state clearly its assumptions. 

We test this method with a range of equations of state, 
including an enhanced Hell photo-heating rate (assumed 
to be due to neglected optical transfer effects) and X- 
ray Compton heating. We show that the method works 
as long as the power spectrum amplitude is sufficiently 
high, so that 034 > 1.6. If the density fluctuations are too 
small, then one important assumption fails: that there be 
a substantial fraction of lines whose width is dominated 
by thermal broadening. When this occurs, it mimics an 
equation of state which is hotter and flatter. 

Very recently, Schaye et al. (1999) independently inves- 
tigated the feasibility of using this method. They used a 
different method for identifying the Nm — b cutoff, but 
came to quite similar conclusions as presented here, with 
one important exception. They argued that there was 
no cosmological dependence on the cutoff in the iVni — b 
plane, in disagreement with the results in this paper. 
However, they examined a relatively small number of sim- 
ulations which all had similar values of 034 , mostly com- 
parable to, or somewhat larger than, the critical value 
listed above. In this case, it would be very difficult to 
notice the effect due to power. 

Applying our results to two quasar lines-of-sight with 
mean redshifts of z = 2.7 and z — 3.4, we derive a 
temperature-density relation from these two systems that 
is similar to those found in at least some of the simula- 
tions presented here. We find a temperature of approx- 
imately 16000 K for gas with the mean density rising to 
about 35000 K for an overdensity of six. The uncertainty 
of these numbers is around 20%, however without any 
more information about the value of 034, they must be 
treated as upper limits to the temperature of the gas. If 
we were to assume that the power spectrum criterion is 
satisfied then this represents fairly hot gas compared to 
traditional models. The additional Hell photo-heating is 
sufficient to produce this much heat; however, by itself 
Compton X-ray heating is not. 

We turn now from the JVhi — b minimum to dn/db, 
the distribution of b parameters. Following similar earlier 



t al. 1997) which includes random noise with a signal- 



to-noise ratio of 60. 



work (Hui & Rutledge 1999), we present a simple linear 
argument which shows that the other important parame- 
ter in determining the b parameter distribution is the am- 
plitude of the power spectrum. We demonstrate that sim- 
ulations reproduce this scaling {b me d ~ T 1 ' 2 a 34 ), and 
show that the shape of the b distribution stays nearly in- 
variant to changes in temperature or the power spectrum 
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amplitude. Its median value can be given as a simple 
function of the gas temperature and 034 (at z ~ 3). If we 
use the temperature derived from the ./Vhi — b method, 
this implies that 034 ~ 1.5 (with an uncertainty of about 
20%), a value which is in reasonable agreement with other 
methods of determining the amplitude of the power spec- 
trum. It should be kept in mind that since the tempera- 
ture measurement is really an upper limit, then this value 
for (T34 is also an upper limit. 

Theuns et al. (1999) also investigated the cosmological 
dependence of the 6-distribution. They found the results 
sensitive to the gas temperature (as we have here), but 
did not find the sensitivity to power spectrum discussed 
here. Again, it seems likely that this is due to the small 
range in spectral amplitudes of their simulations. 

The degeneracy described in this paper between power 
and temperature means that the b distribution alone will 
not be sufficient to determine the equation of state of 
the gas. This is unfortunate because the evolution of the 
temperature-density relationship can provide constraints 
on other cosmological interesting events. For example, 
if the gas were to be colder above z = 3 (as the data 
presented here might be indicating), one possible expla- 



nation would be the late reionization of helium ( Ftcimcrs 



1 1997| ; JHachnclt fc Stcinmctz 1998J ; [Abel fc Hachnelt 



ct 

19E 9|). This degeneracy can be broken by using other as- 



pects of the Lya forest (e.g. Croft et al. 1998; Machacek 



et i,l. 1999) to independently fix the amplitude of fluctu- 



ations at these scales and redshifts. 
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